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ABSTRACT 

We study the cusp/core problem using a secondary infall model (SIM) that 
takes into account the effect of ordered and random angular momentum, dynami- 
cal friction and baryons adiabatic contraction. The model is applied to structures 
on galactic scales (normal and dwarfs spiral galaxies) and on clusters of galaxies 
scales. Our analysis suggest that angular momentum and dynamical friction are 
able, on galactic scales, to overcome the competing effect of adiabatic contrac- 
tion eliminating the cusp. The slope of density profile of inner haloes flattens 
with decreasing halo mass and the profile is well approximated by a Burkert's 
profile. In order to obtain the NFW profile, starting from the profiles obtained 
from our model, the magnitude of angular momentum and dynamical friction 
must be reduced with respect to the values predicted by the model itself. The 
rotation curves of four LSB galaxies from Gentile et al. (2004) are compared to 
the rotation curves obtained by the model in the present paper obtaining a good 
fit to the observational data. The time evolution of the density profile of a galaxy 
of 10*^ — IO^Mq shows that after a transient steepening, due to the adiabatic con- 
traction, the density profile flattens to a ~ 0. On cluster scales we observe a 
similar evolution of the dark matter density profile but in this case the density 
profile slope flattens to a ~ 0.6 for a cluster of ~ IO^^Mq. The total mass profile, 
differently from that of dark matter, shows a central cusp well fitted by a NFW 
model. 

Subject headings: cosmology: theory - large scale structure of universe - galaxies: 
formation 



1. Introduction 

The formation and structure of dark matter (DM) haloes around galaxies and clusters of 
galaxies poses one of the great challenges to theories of structure formation in an expanding 
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universe. The basic problem of the coUisionless collapse of a spherical perturbation in an 
expanding background was first addressed in the two seminal papers by Gunn & Gott (1972) 
and Gunn (1977), where the cosmological expansion and the role of adiabatic invariance were 
first introduced in the context of the formation of individual objects. Self-similar solutions 
were found by Fillmore & Goldreich (1984) (hereafter FG84) and Bertschinger (1985) who 
found a profile of p oc r~^'^^. Hoffman & Shaham (1985) (hereafter HS), showed that 
if structures form from density extrema the density profile is p oc r~°' with a — 
recovering Bertschinger's (1985) profile for n = and Q — 1. Hoffman (1988) refined the 
calculations of HS and made a detailed comparison of analytical predictions of the SIM with 
the simulations of Quinn, Salmon & Zurek (1986) and Quinn & Zurek (1988). 

In the attempt to relax the restrictive assumptions underlying the analytic investiga- 
tions, to incorporate the full range of non-linear gravitational effects, and to study the role 
of initial conditions in shaping the final structure of the dark matter halos, Quinn, Salmon & 
Zurek (1986) pioneered the use of N-body simulations to study halo formation. While Quinn 
et al. (1986) and Efstathiou et al. (1988) found a connection between the density profiles 
of collapsed objects and the initial fluctuation spectrum for Einstcin-de Sitter universes (in 
particular Efstathiou et al. (1988) found density proflles steepening with increasing spectral 
index n), West et al. (1987) arrived at the opposite conclusion. In any case, the previ- 
ous studies showed that the mass density profiles steepen with decreasing in agreement 
with the result of HS. More recent studies (Voglis et al. 1995; Zaroubi, Naim & Hoffman 
1996) showed a correlation between the profiles and the final structures. Finally Dubinski & 
Carlberg (1991), Lemson (1995), Cole & Lacey (1996), Navarro et al. (1996, 1997) (NFW), 
Moore et al. (1998), Jing & Suto (2000), Klypin et al. (2001), Bullock et al. (2001), Power 
et al. (2003) and Navarro et al. (2004) found that although the spherically-averaged density 
profiles of the N-body dark matter halos are similar, regardless of the mass of the halo or 
the cosmological model, their profiles arc significantly different from the single power laws 
predicted by the theoretical studies. The N-body profiles are characterized by an decline 
at large radii and a cuspy profile of the form p(r) oc r~°, where a < 2 near the center. The 
actual value of the inner density slope a is a matter of some controversy, with NFW sug- 
gesting a — 1, but with Moore et al. (1998), Ghigna et al. (2000) and Fukushige & Makino 
(2001) arguing for a = 1.5, while Jing & Suto (2000) and Klypin et al. (2001) claimed that 
the actual value of a may depend on halo mass, merger history, and substructure. Power 
et al. (2003) pointed out that the logarithmic slope becomes increasingly shallow inwards, 
with little sign of approaching an asymptotic value at the resolved radii. In that case, the 
precise value of a, at a given cut-off scale, would not be particularly meaningful. This result 
has been later confirmed by Hayashi et al. (2003) and Fukushige et al. (2004), and it is 
predicted by several analytical models (e.g., Taylor & Navarro (2001), Hoeft et al. 2003). 
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Finally, Navarro et al. (2004) proposed a new fitting formula having a logarithmic slope that 
decreases inward more gradually than the NFW profile. 

The previous different results show that there is no consensus amongst N-body prac- 
titioners on the shape of the density profile. This unsettled state of affairs illustrates the 
difficulties associated with simulating the innermost structure of CDM halos in a rehable 
and reproducible manner. The high density of dark matter in such regions demands large 
numbers of particles and fine time resolution, pushing to the limit even the largest super- 
computers available at present. 

While numerical simulations universally produce a cuspy density profile, observed rota- 
tion curves of dwarf spiral and LSB galaxies seem to indicate that the shape of the density 
profile at small scales is significantly shallower than what is found in numerical simulations 
(Flores & Primak 1994; Moore 1994; Burkert 1995; Kravtsov et al. 1998; Salucci & Burkert 
2000; Borriello & Salucci 2001; de Blok et al. 2001; de Blok & Bosma 2002; Marchesini et 
al. 2002; de Blok 2003; de Blok, Bosma & McGaugh 2003). It seems that the data generally 
favor logarithmic density slopes close to 0.2 (de Blok 2003; de Blok, Bosma & McGaugh 
2003). Using A^-body simulations Kleyna et al. (2003) showed that Ursa Minor dSph would 
survive for less than 1 Gyr if the DM core were cusped. Additionally, Magorrian (2003) 
found a = 0.551^;^ for the Draco dSph. Gentile et al. (2004) (and similarly Gentile et al. 
2006) decomposed the rotational curves of five spiral galaxies into their stellar, gaseous and 
dark matter components and fitted the inferred density distribution with various models and 
found that models with a constant density core arc preferred. The largest part of studies of 
galaxies arrive at the conclusion that cores are preferred to cuspy profiles. Two exceptions 
are the papers of van den Bosch & Swaters (2001) and Swaters et al (2003). van den Bosch 
& Swaters (2001) studied the rotation curves of 19 dwarfs, and claimed that CDM halos are 
consistent with data. However, as pointed out by Moore (2001), to justify this claim "they 
had to throw away half of the Galaxies and adopt unphysicsl (zero) L/M ratios". Spekkens 
et al. (2005) derived inner dark matter halo density profiles for a sample of 165 low-mass 
galaxies using rotation curves obtained from high-quality, long-slit optical spectra assuming 
minimal disks and spherical symmetry. They measure median inner slopes ranging from 
a = 0.22 ± 0.08 to 0.28 ± 0.06 for various subsamples of the data. 

The discrepancy between smulations and observations has been often signaled as a gen- 
uine crisis of the CDM scenario and has become known as the "cusp/core" problem. Since 
LSB galaxies are thought to be ideal for the comparison with theory, as their dynamics are 
dominated by dark matter with little contribution from baryons (Bothun et al. 1997), the 
discrepancy with simulations is particularly troublesome. The significance of this disagree- 
ment, though, remains controversial and different solutions have been proposed. A number 
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of authors attribute the problem to a real failure of the CDM model, or to that of simula- 
tions (de Blok et al. 2001a; de Blok, McGaugh, & Rubin 2001b; Borriello & Salucci 2001; 
de Blok, Bosma, & McGaugh 2003). This has led to suggestions that dark matter proper- 
ties may deviate from standard CDM and several alternatives have been suggested, such as 
warm (Colin et al. 2000; Sommer-Larsen & Dolgov 2001), repulsive (Goodman 2000) fluid 
(Peebles 2000), fuzzy (Hu et al. 2000), decaying (Gen 2001), annihilating (Kaplinghat et 
al. 2000), or self-interacting (Spergel & Steinhardt 2000; Yoshida et al. 2000; Dave et al. 
2001) dark matter. Others argue that the inconsistency may reflect the finite resolution of 
the observations that has not been properly accounted for in the analysis of the HI rotation 
curves (van den Bosch et al. 2000; van den Bosch Sz Swaters 2001). Alternatively, it has 
been suggested that stellar feedback from the first generation of stars formed in galaxies 
was so efficient that the remaining gas was expelled on a timescale comparable to, or less 
than, the local dynamical timescale. The dark matter subsequently adjusted to form an 
approximately constant density core (e.g., Gelato & Sommer-Larsen 1999). This is however 
unlikely to affect cluster cusps. 

On cluster scales. X-ray analyses have led to wide ranging results, from a = 0.6 (Ettori 
et al. 2002) to a = 1.2 (Lewis et al. 2003) or even a = 1.9 (Arabadjis et al. 2002). Mea- 
surements based on gravitational lensing yield conflicting estimates as well, either in rough 
agreement with the results of numerical simulations (e.g. Dahle et al. 2003; Gavazzi 2003), 
or finding much shallower slopes, a — 0.5 (e.g. Sand et al. 2002; Sand et al. 2004). Ricotti's 
(2003) N-body simulations suggest that density profile of DM haloes is not universal (in 
agreement with Jing & Suto 2000; Subramanian et al. 2000; Simon et al. 2003b; Gen et 
al. 2004; Ricotti & Wilkinson 2004; Ricotti et al. 2007), presenting shallower cores in dwarf 
galaxies and steeper cores in clusters. Thus, it seems that there is some evidence coming 
from observations and simulations for a dependence of the inner slope of dark matter haloes 
on halo mass, and hence for the non-universality of the dark matter profile. A conclusive 
theoretical prediction of the central mass distribution of CDM haloes is therefore an impor- 
tant check for any model of structure formation. 

The controversy regarding the "universal" density profile and its logarithmic slope at the 
centre has stimulated a great deal of analytical work. Apart the theoretical work previously 
described, starting fron Gunn & Gott (1972) model, modifications of the self-similar collapse 
model to include more realistic dynamics of the growth process have been proposed (e.g. 
Avila-Reese et al. 1998; Nusser & Sheth 1999; Henriksen & Widrow 1999; Subramanian 
et al. 2000; Del Popolo et al. 2000 (hereafter DP2000)). Several authors (e.g. Syer & 
White 1998; Salvador-Sole et al. 1998; Manrique et al. 2003) argue that the central density 
profile is linked to the merging history of dark matter substructure, and baryons have been 
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invoked both to shallow (El-Zant et al. 2001 (hereafter EZOl), 2004; Romano-Diaz et al. 
2008) and to steepen (Blumenthal et al. 1986) the dark matter profile. Numerous authors 
have emphasized the effect of an isotropic velocity dispersion (thus of non-radial motion) in 
the core of collisionless haloes. Ryden & Gunn (1987) (hereafter RG87) were the first to 
relax the assumption of purely radial self-similar collapse by including non-radial motions 
arising from secondary perturbations in the halo. Gurevich & Zybin (1988a,b) developed 
a formahsm based on the theory of adiabatic capture (e.g., Lifshitz & Pitaevskii 1981), 
to estimate the mass profile in a spherical collapse including non-radial motions. White 
& Zaritsky (1992) introduced a heuristic source term that switches off at turnaround, or 
in another context simply by assigning an angular momentum distribution at turn-around 
time (Sikivie et al 1997). In Avila-Reese et al. (1998) dark matter particles were endowed 
with "thermal motions" resulting in non-radial velocities. Huss et al. (1999b) showed an 
NFW-type density cusp fiattening relative to the isothermal profile just where the velocity 
dispersion changes from predominantly radial to isotropic. This fiattening did not appear, 
in Huss et al (1999a), for the case of pure radial force. A similar result was obtained by 
Tormen et al. (1997), and Teyssier et al. (1997). Nusser (2001) tried two different analytical 
schemes for non-radial velocities in his halos, and concluded that angular momentum is 
most effective when added to a particle at the time of maximum expansion. Hiotelis (2002) 
found that larger amount of angular momentum leads to shallower final density profiles in 
the inner region. EZOl showed that if the gas is distributed in clumps, dynamical friction 
acting on these clumps moving in the background of dark matter particles, dissipate the 
clumps orbital energy and deposit it in the dark matter with the final effect of erasing the 
cusp. Le Delliou & Enriksen (2003) studied the effects of angular momentum on density 
profile in the SIM showing that angular momentum induces an inner turn around radius 
at the size of the self-similar core. Particles with smaller angular momentum will be able 
to enter the core but with a reduced radial velocity compared with the purely radial SIM. 
Ascasibar et al. (2004) included non-radial motions obtaining the result that dark matter 
profile is entirely determined by the initial conditions. Williams et al. (2004), using a 
semi-analytic scheme based on RG87, explored the relationship between the specific angular 
momentum distribution in a halo and its density profile. Compared to those formed in N- 
body simulations, their "semi-analytic" halos are more extended, have fiatter rotation curves 
and have higher specific angular momentum. 

Simulations and semi-analytic models agree on outer parts of the haloes' structure and 
also on properties of haloes substructure; the disagreement in the inner regions could be 
connected to limits in numerical simulations (de Blok 2003; Taylor et al. 2004; see also 
Section 3 of the present paper for a discussion) or to the fact that dissipationless N-body 
simulations does not take into account the effects of baryons on dark matter evolution (see 
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Section 3) . Interestingly, the amount of central substructure seen in the semi-analytic haloes 
is consistent with the amount of substructure inferred from strong lensing experiments (Tay- 
lor et al. 2004). Thus the semi-analytic haloes might provide a more accurate picture of the 
spatial distribution of substructure around galaxies even if an analytical method, no matter 
how sophisticated, will never be able to capture the full extent of complexity of a non-linear 
process. 

Quite apart from the apparent discrepancies between numerical results and observa- 
tions, there is yet another problem, on which very little light has been shed in spite of the 
tremendous progress achieved by numerical work: namely the problem of the physics under- 
lying the universal density profiles predicted by simulations. Intuitively, one could argue that 
the major difference between earlier analytical work and numerical simulations is that halo 
formation in N-body simulations proceeds through repeated mergers. Syer & White (1998) 
and Nusser & Sheth (1999) claimed that the universal profile is a result of hierarchical clus- 
tering by mergers of smaller halos into bigger ones. However, Moore et al. (1999) performed 
N-body simulations with a cut-off in the power spectrum at small scales and also obtained 
halos with cuspy density profiles. This proves that merging and substructure does not play 
a critical role in the formation of density cusps. Huss et al. (1999a,b) found that simulations 
of isolated halos collapsing more or less spherically also result in universal profiles, thus sug- 
gesting that hierarchical merging is not crucial to the outcome. Instead, they suggested that 
the profile is a consequence of a near universal angular momentum distribution of the halos. 
Unfortunately, it is unclear what circumstances lead the halos to this universal angular mo- 
mentum distribution. Thus the issue remains unsolved. Nevertheless the very large amount 
of work carried out by many researchers to date using N-body simulations has met with lim- 
ited success in elucidating the physics of halo formation. The reason is due to the fact that 
the point of force of numerical simulations (namely to capture the full extent of complexity 
of a non-linear process) is also their weakness: numerical simulations yield little physical 
insight beyond empirical findings precisely because they are so rich in dynamical processes, 
which are hard to disentangle and interpret in terms of underlying physics. Analytical and 
semi-analytical models are much more flexible than N-body simulations (see WiUiams et al. 
2004) . So even if analytical models hke SIM treat collapse and viriahzation of halos that are 
spherically simmetric, that have suffered no major mergers, and that have suffered quiescent 
accretion, they are worth investigating (see next Sections for a discussion). 

In this paper, I shall present an analytical model for haloes formation based on SIM. 
The paper is an extension of DP2000 in which we showed that a simple spherical infall model 

gives rise to profiles which are not power laws. So, again, in the present paper, we will follow 
the HS spirit assuming that objects form around maxima of the smoothed density field and 
the model of Zaroubi and Hoffman (1993) (hereafter ZH93). Differently from DP200 and 
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ZH93, we shall study the collapse in presence of non-radial motions (ordered and random 
angular momentum), dynamical friction and baryons adiabatic contraction. 

The plan of the paper is the following: in Section 2, we introduce the model that shall 
be used to calculate the density profiles. Section 3 deals with results and discussion. Section 
4 is devoted to conclusions. Appendices describe the initial set-up, namely the way initial 
conditions, angular momentum, dynamical friction and adiabatic contraction are introduced 
in the model and, how the adiabatic contraction is calculated. 

2. Model 

The simplest version of SIM considers an initial point mass, which acts as a nonlinear 
seed, surrounded by a homogeneous uniformly expanding universe. Matter around the seed 
slows down due to its gravitational attraction, and eventually falls back in concentric spher- 
ical shells with pure radial motions. The assumptions of SIM that are most often questioned 
are the spherical symmetry and the absence of peculiar velocities (non-radial motions): in 
the "real" collapse, accretion does not happen in spherical shells but by aggregation of 
subclumps of matter which have already collapsed; a large fraction of observed clusters of 
galaxies exhibit significant substructure (Kriessler et al. 1995). Motions are not purely ra- 
dial, especially when the perturbation detaches from the general expansion. Nevertheless 
the SIM gives good results in describing the formation of dark matter haloes, because in 
energy space the collapse is ordered and gentle, differently from the chaotic collapse that is 
seen in N-body simulations (Zaroubi, Naim & Hoffman 1996). This is confirmed in other 
studies (Toth & Ostriker 1992; Huss, Jain & Steinmetz 1999a,b; Moore et al. 1999). More- 
over, judging by the commonness of extended thin spiral disks in the Universe, Dark halo 
formation of disk galaxies may reasonably be described by the SIM because dynamical disk 
fragilty implies that major mergers could not have played a significant role in these cases 
(e.g., Toth & Ostriker 1992). We should also add that analytical and semi-analytical models 
have some advantages on N-body simulations: a) they are flexible (one can study the effects 
of physical processes one at a time); b) one can incorporate many physical effects at least in 
a schematic manner; c) they are computationally efficient (it takes about 10 s to compute 
the density profile of a given object at a given epoch on a desktop PC (Ascasibar et al. 
2007). 

As I showed in DP2000, the discrepancies between the SIM and some high resolution 
N-body simulations are not due to the spherical symmetry assumption of the SIM but arises 
because of some non-accurate assumptions used in its implementation. In the quoted paper, 
I improved HS model showing that SIM predicts non power-law profiles. HS considered 
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a scale-free initial perturbation spectra, P{k) oc fc" (where n is the spectral index) and 
assumed that local density extrema are the progenitors of cosmic structures and that the 
density contrast profile around maxima is proportional to the two-point correlation function. 
They thus showed that p oc with a = 3(3 -|- ?t,)/(4 -|- ra), being n the spectral index. 

The modification of HS model in DP2000 are connected to the following observations: 
a) the conclusion that the final density of the profile is p oc for n < —1, claimed by HS, 
is not a direct consequence of the SIM model, but it is an assumption made by the quoted 
authors, following the study of self-similar gravitational collapse by FG84. b) The assumption 
made by HS that the initial mean fractional density excess inside the shell, 6{r) = 
is proportional to the correlation function f(r) (namely 5{r) oc ^(r) oc r~^^~^"'^) is not good 
for regions internal to the virial radius, |3 (see Peebles 1974; Peebles & Groth 1976; Davis 
& Peebles 1977; Bonometto & Lucchin 1978; Peebles 1980; Fry 1984; DP2000). In other 
words, HS's solution applies only to the outer regions of collapsed haloes, and consequently 
the conclusion, obtained from that model, that dark matter haloes density profiles can be 
approximated by power-laws on their overall radius range is not correct, c) According to 
Bardeen et al. (1986), (hereafter BBKS), the mean peak profile depends on a sum involving 
the initial correlation function, ^(r) oc r' 

-{3+n)^ and its Laplacian, V^^(^) r-(5+") (BBKS; 
RG87) (see Appendix B), not only the correlation function ^(r) as in HS. 

As shown in DP2000 (their Eq. 20, and Fig. 6 of the present paper), the initial mean 
density obtained using the model of that paper is extremely different from that obtained 
and used in HS, and differently from HF's model the density profiles are not power-laws but 
have a logarithmic slope that increase from the inner halo to its outer parts. 

In summary the model assumes that the initial probability distribution of the density 
field is Gaussian. The dynamical evolution of matter at the distance Xi from the peak is 
determined by the mean cumulative density perturbation within Xi (see Appendix A Eq. 
flA2p ) and the maximum radius of expansion can be obtained knowing Xj and the mean 
cumulative density of the perturbation (Eq. flA2p ). After reaching maximum radius, a shell 
collapses and will start oscillating and it will contribute to the inner shells with the result 
that energy will not be an integral of motion any longer. The dynamics of the infalling shells 
is obtained by assuming that the potential well near the center varies adiabatically (Gunn 
1977, FG84). The details of the model, which is an extension of DP2000 to take account 



'^Pb is the mean background (critical) density. 

^The virial radius of the halo is defined by the radius of a sphere enclosing a given over-density, Ay. This 
last value is obtained in our case from Bryan & Norman (1998). The relation between the virial radius and 
the virial mass is given by = 47r/3/9bA^rj] 
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of the effects of angular momentum, dynamical friction and dissipative infall on the density 
profile, are exposed in Appendix A. 

The quoted model needs that the initial density profile is given together with angular 
momentum and dynamical friction. The way these quantities are calculated and introduced 
in the model are described in Appendix B, C and D, respectively. Appendix E shows how 
we take into account adiabatic contraction. 



3. Results and discussion 

After fixing the initial conditions and describing how to calculate angular momentum, 
dynamical friction and adiabatic contraction, we can use the model in Section 2 to obtain 
the density profile of haloes. As reported in Appendix E, we employ the usual 
assumption that initially baryons had the same density profile as the dark matter 
(Mo et al. 1998; Cardone & Sereno 2005; Treu & Koopmans 2002; Keeton 2001; 
Klypin 2002; Tonini et al. 2006). In panels (a)-(d) of Fig. 1, the solid line represents 
the NFW profile for haloes having masses equal to lO^M© (panel a), lO^^M© (panel b), 
10^-^ Mq (panel c), IO^^Mq (panel d). The NFW profile for the given mass was calculated by 
means of the relationships connecting the concentration parameter, c, and the virial mass. 
My, to the shape of NFW profile. We used the following equation for c: 

(Gentile et al. 2007), and the usual one for the NFW profile 

o(r) = = (2) 

^ r/r,(l + r/r,)2 r/r,(l + r/r,)^ 

where 

= "#log(l + c)-c/(l + c) 

and A^, is the virial overdensity (see Bryan & Norman 1998). The scahng radius, r^, of NFW 
profile is connected to the virial radius, concentration parameter and virial mass through 
c = Ty/rg, where: 

(Gentile et al. 2007). 
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Notice that we took the redshift dependence in the model using the technique 
described in Del Popolo (2001) (the reader is referred to the quoted paper to 
have more insights). 

NFW profiles change slope rapidly from a = 1 to a = 3 at the characteristic radius r^. 
For example, in the case of the 10^"^ Mq, the NFW profile have a characteristic scale length, 
equal to 0.1r„, beyond which the density profile steepens, so that much of the mass is piled up 
within 10% of the virial radius. The haloes obtained using the model in Section 2 are different 
in character from the profiles predicted by numerical simulations, like those of NFW. Within 
the virial radius the log-log density slope changes gradually and the slopes of the inner part 
of haloes flattens with decreasing mass. The dotted line represents the density profile of 
haloes calculated according to the model of Section 2 and have masses: 10^ Mq (panel a), 
1O^°M0 (panel b), lO^^M© (panel c), lO^^M© (panel d). Our results shows a steepening of 
the density profile with increasing mass with slopes a ^ for M ^ 10^ — 10^ Mq, a ~ 0.2 for 
M ~ lO^^M©, a ~ 0.6 for M ~ IO^Mq, a ~ 0.8 for M ^ W^^Mq. The dashed line in Fig. 
1 (almost indistinguishable from the density profile of the IO^Mq halo), represents a fit to 
the density profile by means of a Burkert's profile considered a good fit to the dark matter 
rotation curves inferred from observations (e.g., Salucci & Burkert 2000). The functional 
form of this profile is characterized by: 

^^'^^ = (l + r/r,)[l + (r/r„)2] 

where po — ps and — (EZOl). The dark matter within the core is given by Mo = 
LGpoTo- Although the dark matter parameters Tq, po and Mo are in principle independent, 
the observations reveal a clear connection (Burkert 1995): 

Mo = 4.3 X 10^ {ro/kpcf^ Mq (6) 

which indicates that dark haloes represent a one parameter family that is completely speci- 
fied, e.g., by the core mass (Salucci & Burkert 2000). 

Two important things must be noticed: a) less massive haloes are less concentrated; b) 
the halo's inner slope is smaller for smaller mass. 



The first point can be explained as follows: higher peaks (larger z/), which are progeni- 
tors of more massive haloed, have greater density contrast at their center, and so shells do 



^ This affirmation is discussed in Peacock &c Heavens 1990; Del Popolo & Gambera 1996 
or Gao &; White (2007) (Fig.l). Modelling the peaks as triaxial spheroids one obtains a peak 
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not expand far before beginning to collapse. This reduces j and allows haloes to become more 
concentrated. An alternative explanation is connected to the quoted angular momentum- 
density anti-correlation showed by Hoffman (1986): j oc z/~^/^. So, density peaks having 
low (high) value of u acquire a larger (smaller) angular momentum than high u peaks and 
consequently the halo will be less (more) concentrated. It is important to notice that the 
quoted trend of increased central concentration as a function of mass applies only to halos 
that started out as peaks in the density field smoothed with a fixed Rf scale. Our conclu- 
sions do not mean that, for example, clusters of galaxies will be very much more centrally 
concentrated than galaxies, since different smoothing scales would apply in the two cases. 

Point (b) can be explained in a similar way to (a), as described previously. Less massive 
objects are generated by peaks with smaller u, which acquire more angular momentum [h 
and j). Angular momentum sets the shape of the density profile at the inner regions. For 
pure radial orbits, the core is dominated by particles from the outer shells. As the angular 
momentum increases, these particles remains closer to the maximum radius, resulting in a 
shallower density profile. Particles with smaller angular momentum will be able to enter 
the core but with a reduced radial velocity compared with the purely radial radial SIM. For 
some particles the angular momentum is so large that they will never fall into the core (their 
rotational kinetic energy makes them unbound). Summarizing, particles with larger angular 
momenta are prevented from coming close to the halo's center and so contributing to the 
central density. This has the effect of flattening the density profile. This result is in agreement 
with the previrialization conjecture (Peebles & Groth 1976; Davis & Peebles 1977; Peebles 
1990), according to which initial asphericities and tidal interactions between neighboring 
density fluctuations induce significant non-radial motions that oppose the collapse. In order 
to reproduce the NFW profile, we performed an experiment similar to that performed by 
Williams et al. (2004), namely we reduced the magnitude of the h and j angular momentum 
and dynamical friction, /x. The experiment was performed on the halo of mass lO^^M©, 
and in order to reproduce the NFW profile having c = 10 and mass ~ lO^^M©, we had 
to reduce the magnitude of h of a factor of 2, j and yU of a factor 2.5. The result of the 
quoted experiment is the dashed line in Fig. Id, which closely reproduces the NFW profile. 
Similarly, Williams et al. (2004) had to reduce random velocities, which amount to reducing 
the angular momentum, in order to obtain a NFW profile. With each reduction of the 
random velocities, the profiles get steeper at the center. This effect can be understood, as 



mass M = '^^a_l:'^Qgf^^^''3%' for 0.5 < 7 < 0.8, where i?* and 7 are given in Eq. (|B6|) and Eq. (fBT]) . 
showing an increase of mass with ly. It is reasonable that lower v peaks should have lower 
mass; peaks with ~ will tend to sit in regions of larger scale underdensity (cancelling the 
small-scale overdensity), and hence the ~ PbRl of material, which initially surrounds the peak, 
may not be accreted following central collapse. 
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already reported, as follows: the central density is built up by shells whose pericenters are 
very close to the center of the halo. Particles with larger angular momenta are prevented from 
coming close to the halo's center and so contributing to the central density. The correlation 
between increasing angular momentum and the reduction of inner slopes in halos has been 
also noticed by several other authors (Avila- Reese et al. 1998, 2001; Subramanian et al. 
1999; Nusser 2001; Hiotelis 2002; Le Delhou & Henriksen 2003; Ascasibar et al. 2003). 

Before going on, we want to add that haloes of a given mass corresponding 
to higher peaks (larger u) in the initial density profile, collapse at earlier epochs 
than those having smaller u and they give rise to more massive and concentrated 
haloes. Higher peaks, being more concentrated, feel less the tidal torque and 
as a consequence the central part of the profile is steeper than in peciks having 
smaller u which are more torqued and collapse in flatter profiles. 

As Shown in Fig. 2, the solid histogram representing the total specific an- 
gular momentum distribution of the density profile reproducing the NFW halo 
(described in the previously quoted experiment) is more centrally concentrated 
than the total specific anguleir momentum distribution of our reference haloes 
(dashed histogram), and is closer to those of typical halos emerging from nu- 
merical simulations. In Fig. 2, the dotted-dashed and dashed line represents 
the quoted distribution for the halo n. 170 and n. 081, respectively, of van den 
Bosch et al. (2002). The halo n. 170 resembles most of the specific angular mo- 
mentum distributions, while the halo n. 081 has the shallowest distribution in 
their simulations. This may suggests, in agreement with Williams et al. (2004), 
that haloes in N-body simulations lose a considerable amount of anguleir mo- 
mentum between 0.1 and 1 r^. Since virialization proceeds from inside out, this 
means that the angulEir momentum loss teikes place during the later stages of the 
halos' evolution, rather then very eeirly on. This is somehow confirmed by the so 
called angular momentum catastrophe, namely the fact that dark matter halos 
generated through gas-dynamical simulations are too small and have too little 
angular momentum compared to the halos of real disk galaxies, possibly because 
it was lost during repeated collisions through dynamical friction or other mech- 
anisms (van den Bosch et al. 2002; Naveirro &; Steinmetz 2000). The problem 
can be solved invoking stelleir feedback processes (Weil et al. 1998), but part 
of the anguleir momentum problem seems due to numerical effects, most likely 
related to the shock capturing, cirtificial viscosity used in smoothed peirticle hy- 
drodynamics (SPH) simulations (Sommer-Leirsen &;Dolgov 2001). We discussed 
the effect of changing the magnitude of angular momentum but wc did not speak of the 
effect of changing the magnitude of /i (dynamical friction) . The effect of changing this last 
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quantity is very similar to changing the magnitude of angular momentum: an increase in 
the term fi produces shallower profiles as larger values of angular momentum does. This is 
expected from Fig. 11, showing that dynamical friction influence the dynamics of collapse 
in a similar way to that of angular momentum slowing down the collapse of outer shells and 
so compelling the particles to remain closer to the maximum radius. 

Fig. 3 shows the evolution of a IO^Mq halo. The evolution was obtained by 
calculating the profile at different redshifts (see Del Popolo 2001). We started 
the evolution at z = 50. The solid line represents the profile at z = 10. This is the 
epoch at which the profile virialize^ and the shape of the profile at virialization 
is obviously given by the line corresponding to z = 10. At subsequent times, the 
profile ai z = 3, z = 2, z = 1 and z = is represented by the long-dashed line, short-dashed 
line, dot-dashed line and dotted-line, respectively. The evolution after virialization 
is produced by secondary infall, two-body relaxation, dynamical friction and 
angular momentum. The cusp is slowly eliminated and within ~ 1 kpc a core forms. At 
an early redshift, 2 ~ 5, the dark matter density experiences the adiabatic contraction by 
baryons producing a slightly more cuspy profile than that represented by the solid line (not 
represented in the plot). 

The previous result is similar to what found by Romano-Diaz et al. (2008) who studied 
the dark matter cusp evolution using N-body simulations with and without baryons. The 
"erasing" of the cusp is associated by them to the heating up of the cusp region via dynam- 
ical friction (EZOl) or influx of subhaloes into the innermost region of the dark matter halo. 
In our model the erasing of the cusp is connected to the joint effect of dynamical friction and 
angular momentum. As previously discussed, larger amount of angular momen- 
tum leads to shallower final density profiles in the inner region because particle 
with larger angular momentum have lower probabilities to enter the center. The 
effects of dynamical friction can be interpreted in two different fashions: (a) 
an increase in the term fi is very similar to changing the magnitude of angular 
momentum (see Fig. 11) with the final result of producing shallower profiles; 
(b) dynamical friction can act on gas moving in the background of dark matter 
particles, dissipate the clumps orbital energy and deposit it in the dark matter 
with the final effect of erasing the cusp (similarly to EZOl; El-Zant et al. 2004; 
TLS; Romano-Diaz et al. 2008;). Baryons have another effect, at an early red- 
shift, the dark matter density experiences the adiabatic contraction by baryons 
producing a slightly more cuspy profile. This last is overcome from the previous 



The virialization is calculated as usual, by means of the virial theorem, as the epoch at which the ratio 
between the kinetic and potential enrgy is ~ 1/2. 
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two effects. As shown by Fig. 11, the magnitude of dynamical friction effect is 
a bit larger than that due to angular momentum and that these two effects add 
to improve the flattening of the profile. 

In Fig. 4, we plot the rotation curves obtained by our model and we compare them to 
four LSB galaxies studied by Gentile et al. (2004), namely ESQ 116-G12, ESQ 79-G14, ESQ 
287-G13 and NGC 1090. The physical parameters of the quoted galaxies are given in Gentile 
et al. (2004) (Table 1). In the case of ESQ 116-G12 the total HI mass is Mhi ^ 1.5 x lO^M© 
and the dynamical mas^ M^yn — 3.3 x IO^^Mq. For the other three galaxies ESO 79-G14, 
ESQ 287-G13 and NGC 1090, the masses are: Mhi ^ 3.5 x IO^Mq, Mdyn ^ 1.3 x W^Mq; 
Mhi ^ 1.1 X lO^^Mo, Mdyn ^ 1.9 x IO^^Mq; Mhi ^ 8.5 x IO^Mq, Mdyn ^ 1.8 x W^Mq, 
respectively. In all the four cases, the data are compared with the rotation curve obtained 
using our model (solid line) and with rotation curves obtained from NEW profile (dotted 
lines), given by: 

r Hi + cx)-cx/{i + cx) 1 

'^^'^^"'^''1 x[ln(l + c)-c/(l + c)] / 

where x = r/r^ and Vy is the virial velocity!^. Fig. 4 shows that NEW haloes are higher 
than the rotation curves obtained using our model, in which more massive haloes tend to be 
more centrally concentrated and have flatter rotation curves. Less massive haloes are less 
concentrated, ad have slowly rising rotation curves. In contrast NEW rotation curves rise 
very steeply and as a consequence NEW fits to dwarf galaxy rotation curves have too low 
concentration parameters (van den Bosch & Swaters 2001) compared to N-body predictions. 
NEW fails to reproduce velocities and shape of the observed rotation curves, especially in the 
case of ESO 116-G12, ESO 79-G14. NEW haloes predict too high velocities in the central 
part of haloes, and even leaving c as free parameter, instead of using Eq. ([1]), there is no 
appreciable improvement in the fit. Using Eq. ([1]), one obtains very low values of c. The 
result is similar to that described by Gentile et al. (2004): data are much better described 
by core-like profiles, like the Burkert profile generating flatter rotation curves: 



y(r)=('^^)V'^ln 



(1 + r/ro) V (1 + (r/ ro)^) - arctan(r/ r^) > 



Our rotation curves are very similar to those generated by the Burkert profile and the 
residuals and discrepant points for our rotation curves are close to that given in Gentile et 
al. (2004) for the Burkert's fit to their data. 



^The dynamical mass Mdyn is determined at the farthermost radius with data. 

^The value of the characteristic velocity Vu of the halo is defined in the same way as the virial radius r„. 
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As noticed in the introduction, while on galactic scales a large number of studies predicts 
central cores and it seems that the cusp / core problem is a real problem not attributable to 
systematic errors in the data (de Blok, Bosma & McGaugh 2003), on cluster scales the 
situation is less clear with slopes ranging from a — 0.5 (e.g. Sand et al. 2002; Sand et al. 
2004) to a = 1.9 (Arabadjis et al. 2002). 

In order to study the problem on cluster scales, wc have calculated the density profile 
evolution of dark matter and that of the total matter distribution for halo of ~ 10^^h~^MQ. 

Fig. 5 plots the evolution of a density profile of lO^^/i^^M©. The solid line repre- 
sents the density profile at z = 3 which slightly steepens due to baryon settling 
in virialized dctrk matter haloes (AC) at z = 2 (not shown in the plot). Once the 
baryons condense to form stars and galaxies, they experience a dynamical fric- 
tion force from the less massive dark matter particles as they move through the 
halo. Energy and angular momentum is transferred to dark matter, increasing 
its random motion. Moreover, angulcir momentum acquired in the expansion 
phase gives rise to non-radial motions in the collapse phase. The effect of angu- 
lcir momenta and dynamical friction overcomes that of the AC and the profile 
stEirts to fiatten (z — 1.5 dotted line; z — 1 short-deished line; z — long-dashed 
line). The final dark matter profile (long-dashed line) is characterized by a log- 
log slope of a ~ 0.6 at O-Olr^,. So the situation is similar to that of haloes on 
galactic scales but the slope is larger than for dwarf galaxies. In this case the 
profile virializes at 2; ~ and the corresponding line respresents the profile at 
virialization. Comparing the two profiles at virialization in Fig. 1 and Fig. 3, 
(line corresponding to 2; = 10 for the profile having lO^M©; line corresponding to 
2; ~ for the profile having IO^^Mq), it is cleeir that the profile at virialization is 
fiatter in the case of smaller masses. The difference is due to the fact that anguleir 
momentum starts to change the profile shape steirting already from turn-around 
epoch. The dot-dashed line in Fig. 5 (uppermost line), represents the final total density 
profile of a IO^'^Mq halo. The plot shows that the cusp in the total density profile is not 
"erased" as in the case of the dark matter profile. This result also implies that the baryonic 
component becomes steeper than the original NFW profile. The behavior of total mass is 
in agreement with X-ray observations by Chandra and XMM (Buote 2003, 2004; Lewis et 
al. 2003) weak lensing (Dahle et al. 2003) and strong lensing (Bartelmann 2002), which are 
consistent with a cusp having a = 1 or larger. The behavior of the dark matter halo is in 
agreement with analysis of Sand et al. (2002, 2004) who fitted the baryonic and dark matter 
profiles only in the very inner part of the cluster MS 2137-23 within ~ 50h~^ kpc by means 
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of a generalized NFW profile: 



o(r) = (9) 



obtaining a nearly fiat core a = 0.35. The steepening of the baryonic component, is consistent 
with what found by Brunzendorf & Meusinger (1999), who found that the projected galaxy 
distribution in Perseus cluster diverge as r~^. 

The results previously reported have several implications on the effort to test predictions 
of the CDM model observationally. The test that received much attention in the last decade, 
as several times stressed, is the density distribution in the inner regions of galaxies and 
clusters. 

As previously stressed, in our model the cusp/core problem is solved by 
the effects of dynamical friction and angular momentum overcoming that of the 
AC, however, at the same time, the effect of baryons is non-negligible (e.g. when 
compcired to collisionless N-body simulations), in the sense that we will describe. 
The effect of bairyons is two-fold: a) in the inner pcirts of the haloes beiryons Eire 
more present and their adiabatic collapse, from one side, steepens the density 
profile because the amount of dark matter in the central region will always be 
increased, b) From the other side, one has to take account of the exchange of 
angular momentum between the baryons and the dark matter. Dynamical fric- 
tion can result in a transfer of angular momentum from the baryons to the dark 
matter. Because the dark matter gains angular momentum, it moves further 
from the galactic center reducing the steepeness of the inner profile (see also 
Klypin et al. 2001). So beiryons have orbital energy that is transferred through 
dynamical friction and deposited in the deirk matter, giving a pEirtial solution to 
the cusp/core problem (EZOl; El-Zant et al. 2004; Romano-Diaz et al. 2008). 
In collisionless N-body simulations, this complicated interplay between different 
effects is not taken into account and makes it necessary to run N-body simula- 
tions that repeat the mass modeling including a self-consistent treatment of the 
baryons and dark matter component. So, on galactic scales, where dark matter 
dynamics and baryons dynamics are entangled, the cusp/core problem seems to 
be a "genuine" one, in the sense that the disagreement between observations 
and N-body simulations is not due to numerical eirtifacts or problems with sim- 
ulations. At the same time it is an appeirent problem, since the disagreement 
between observations and dissipationless simulations is related to the the fact 
that this last are not taking account of baryons physics. This means that we are 
compciring two different systems, one dissipationless (i.e., DM) and the other 
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dissipational (i.e., inner pcirt of structures), and we cannot expect them to have 
the same behavior. 

As our results show, baryons at early times steepens the cusp due to the AC and after 
the cusp is erased through dynamical friction and non-radial motions effects. This results are 
in agreement with other analyses which studies separately the effects of dynamical friction 

(e.g., EZOl) and those of angular momentum (non-radial motions) (e.g., Nusser 2001; Hiotelis 
2002; Ascasibar et al. 2003; Williams et al. 2004) and with the recent Sph simulations of 
Romano-Diaz et al. (2008). Going to larger scales the situation changes. The analysis of 
density distribution for bright galaxies is complicated by the uncertain contribution of stars 
to the total mass profile (Treu & Koopmans 2002; Mamon & Lokas 2004). Some analyses 
tend to favor inner slopes shallower than predicted by CDM (e.g.. Gentile et al. 2004) 
but other deduce slopes of the inner profiles that are at least marginally consistent with 
predictions (Treu & Koopmans 2002, 204; Koopmans & Treu 2003; Jimenez et al. 2003). 
As previously reported, our results shows a steepening of the density profile with increasing 
mass with a density profile of haloes of mass > lO^^M© having slopes > 0.8. This is in 
agreement with recent N-body simulations having a logarithmic slope that decreases inward 
more gradually than the NFW profile (Hayashi et al. 2003; Navarro et al. 2004; Stadel et 
al. 2008). In the case of Stadel et al. (2008) the logarithmic slope is 0.8 at 0.05% of Vy. 

The density distribution in clusters of galaxies can, in principle, provide a cleaner test 
of the models because the effects of the baryons and gas on the dark matter distribution are 
expected to be smaller and simpler. However also in this case observations predict slopes 
ranging from a < 0.5 (a = 0.35, Sand et al. 2002, 2004) to values larger than one (Arabadjis 
et al. 2002) and in some cases different results even for the same object. It is the case of 
the cluster MS2137-23 studied by Sand et al. (2004), who found a shallow density slope 
(< 0.5) while Dalai & Keeton (2003), Bartelman & Meneghetti (2003) and Gavazzi (2003) 
contested Sand's results, which according to them is neglecting lens ellipticity, and found 
consistency of the inner slope with a NFW profile. Our result concerning cluster scales makes 
a difference between dark matter and total mass distribution: the first tend to be less cuspy 
than observations in agreement with some observations (e.g.. Sand et al. 2002, 2004), while 
the second is cuspy and well described by a NFW profile (in agreement with Brunzendorf & 
Meusinger 1999) . 

One question that may arise at this point is: why the results of analytical or semi- 
analytical models are different from those of N-body simulations? 

The discrepancy between N-body simulations and observations at small scales has led 
some authors to attribute the problem to a real failure of the CDM model, or to that of 
simulations (de Blok et al. 2001a; de Blok, McGaugh, & Rubin 2001b; Borriello & Salucci 
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2001; de Blok, Bosma, & McGaugh 2003; de Blok 203; Taylor 2004). One of the reasons 
why N-body simulations could give unreliable results are connected to two-body relaxation. 
The processes of relaxation is difficult to quantify, but in the large N limit one expects 
that the discreteness effects inherent to the N-body technique vanish, so one tries to use as 
large a number of particles as computationally possible. Unfortunately in most cosmological 
simulations the importance of two body interactions does not vanish if one increases A^, 
since structure formation in the cold dark matter (CDM) model occurs hierarchically since 
there is power on all scales, so the first objects that form in a simulation always contain only 
a few particles (Moore et al. 2001), (Binney & Knebe 2002). With higher resolution the 
first structures form earlier and have higher physical densities because they condense out 
of a denser environment. Two body relaxation increases with density, so it is not clear if 
increasing the resolution can diminish the overall amount of two body relaxation in a CDM 
simulation, i.e. if testing for convergence by increasing the mass resolution is appropriate. 
Diemand et al. (2004a,b) explored the effect of resolution on the degree of relaxation finding 
that increasing N slowly reduces the degree of relaxation oc N'^'"^^ rather than proportional 
to N as expected from the coUisionless Boltzmann equation. This means that if to resolve 
10% of the virial radius we need 1000 particles, we need 10^ particles to resolve 1% of the 
virial radius. To have such a high number of particles per halo the use of the multi-mass 
technique is required. This technique consists of simulating with high mass resolution only 
the particles that will end up in the halo of interest at ^ = 0, while having less mass resolution 
for the particles that end up far away from the halo of interest. The task of achieving reliable 
profiles with such a high spatial resolution is sensitive to numerical integration errors and 
requires careful resolution studies (Power et al. 2003); perhaps this is a reason for the 
disagreement between groups, using different codes, on the slope of the inner profile (Ricotti 
2003). Moreover, this has the drawback that each simulation can resolve only one halo at 
a time, therefore selection effect biases, (difficult to control, determined by the criteria for 
picking the halos to re-simulate) and a poor statistical sample could affect the reliability of 
the final result even if the simulation is very accurate and has high resolution (Ricotti 2003). 
More recently Stadcl ct al. (2008) by means of simulations using several billions particle 
measured the density profile to a distance of 120 pc (0.05% of r^). Convergence in the 
density profile and the halo shape scales as A^~^/^, but the shape converges at a radius three 
times larger at which point the halo becomes more spherical due to numerical resolution. 
This last simulation has surely enough resolution (< 1 kpc) to distinguish between the core 
and cusp model, which was one of the problems de Blok (2003) enumerated in the reasons 
why simulations produce just cuspy profiles. Other problems due to numerical artifacts are 
the over-merging problem, the artificial disruption of substructure due to numerical effects, 
that according to Taylor et al. (2004) has not been solved as claimed by some authors (e.g. 
Ghigna et al. 2000). Even if the effects of particle discreteness in N-body simulations of 
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ACDM are still an intensively debated issue (Romeo et al. 2008), it is highly probable that 
N-body simulations are correctly predicting the density profiles of the CDM haloes but on 
small scales other effects not taken into account by dissipationless simulations, for example 
the presence of baryons and the effect that they have on dark matter, change the halo 
shape. In other words, at this state of the art of N-body simulations and observations, the 
discrepancy between simulations and observations is not the fault of numerical artifacts of 
simulations or problems of the observations (see de Blok (2003) to have a list of the problems 
imputed to observations and its confutation) but it is just connected to the fact that we are 
trying to force dissipationless simulations to predict the same behavior for the density profile 
of a system whose physics is not just the dissipationless physics typic of dark matter. It is 
interesting to note that the problems of CDM only become clear on length scales where the 
baryons start playing a role and that this applies not only to the cusp/core problem but 
also to the missing dwarfs problem. Moreover, it is noteworthy that if, indeed, most star- 
forming galaxies in the early universe lost their DM cusps because of stellar feedback, the 
missing dwarfs (satellites) problem could also be solved. Dwarf galaxies without a central 
cusp have a lower average core density than cuspy ones, and are hence much easier to disrupt 
tidally during the hierarchical assembly of larger galaxies (Mashchenko & Sills 2005). As a 
consequence, the removal of galactic cusps by stellar feedback in the early universe would 
result in fewer satellites today. This again indicates that baryon physics is one of the missing 
pieces of the puzzle, and will very likely make a major contribution toward a solution. If this 
is true, it would be unwise to ignore the conclusions to which data are leading us, namely 
that small scales tells us more about galaxy formation than it doe. about CdB Iu other 
terms, the centers of galaxies are, special places, the only places where we can study dark 
matter under peculiar conditions. 

The predictions of large amounts of small-scale structure and substructure in CDM 
cosmologies is perhaps startling, but it is not in and of itself a reason to reject CDM. 

Another important point on which our results can tell something is the debate on the 
universality of the density profiles of dark matter haloes. The result that the density profiles 
of haloes in CDM and other hierarchical clustering cosmologies have a universal form which 
is well represented by the simple fitting formula given by NFW96, 97, has been confirmed 
in almost all the subsequent papers dealing with the subject except some of them (Jing & 
Suto 2000; Subramanian et al. 2000; Ricotti 2003, Ricotti & Wilkinson 2004; Cen et al. 
2004; Ricotti et al. 2007). Ricotti's papers, using the same type of simulations in NFW97, 
found that at virialization the central logarithmic slopes a at 5%-10% of the virial radius 



Other possibilities are that the observed dark matter cores are teUing us that dark matter has pressure 
at small scales or something unexpected. 
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are correlated with the halo mass, with a = 0.2 for ~ lO^M©, a = 1 for ~ lO^^Mo, a ~ 1.3 
for ^ IO^^Mq, and that there is no reason to believe that the value of a converges to any 
asymptotic value as also suggested by some high-resolution N-body simulations (e.g., Navarro 
et at al. 2004; Graham et al. 2005; Stadel et al. 2008). This leads to the conclusion that 
density profiles do not have a universal shape. Moreover in agreement with Subramanian 
(2000), the halo shape at a given mass or spatial scale depends on the slope of the power 
spectrum at that scale. Similar correlations have been found by other authors (e.g., Jing & 
Suto 2000; Taylor & Navarro 2001; Cen et al. 2004). Cen et al. (2004) confirmed Ricotti's 
result; in addition they identify a redshift dependence of the typical halo profile. Graham et 
al. (2005) and Merrit et al. (2005) also find a correlation between halo mass and the shape 
of the density profile, parameterizing it in terms of the Sersic profile index. There are also 
observational evidences of a mass dependence of the dark matter density profile, as reported 
in Ricotti & Wilkinson (2004) and in the introduction of the present paper. Other examples 
of a possible non-universality of density profiles comes from galaxies observations. Galaxies 
like NGC2976, NGC6689, NGC5949, NGC4605, NGC5963 have very different values of the 
slope: a ~ 0.01, 0.80, 0.88, 0.88, 1.28, respectively (see Simon 2003a,b). Moreover observed 
slopes on galactic scales have large scatter compared to simulations and mean slope shallower 
than simulations. Thus there is also some observational evidence for a dependence on halo 
mass of the inner slope of dark matter haloes, and hence for the non-universality of the dark 
matter profile. 

Our results in Fig. 1 show clearly different slopes flattening for decreasing values of the 
halo mass in agreement with Ricotti (2003), Ricotti & Wilkinson (2004) and Williams et al. 
(2004) (eg., their Fig. 1), showing a flattening of the profile with decreasing mass. As shown 
in Fig. 3, the density profile of a halo of lO^M© at virialization is different from 
a NFW profile, the inner slope is less steep than the NFW model. Dynamical 
friction and angular momentum contribute to flatten the density profile which 
is well fitted by a Burkert's profile at 2; = 0. In the case of the IO^^Mq halo, as 
shown in Fig. 5, the situation is similar except that the final density profile is 
more steep than in the case of a IO^Mq halo. 

The idea coming out from the previous arguments is that one should not 
expect universality of density profile^"]. Summarizing, the flattening of the in- 
ner slopes of haloes is produced by the role of angular momentum, dynamical 
friction and the interplay between dark matter and baryonic component. We 
want to remark that the Aquarius Project showed that even in N-body simulations 



^^As previously stressed, the situation is different for the total mass profile for whieh the NFW seems to 
give a good fit. 
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of ACDM haloes the mass profiles of haloes is not strictly universal (Navarro et 
al. 2008). Concluding, concerning the cusp/core problem, there is no real contradiction 

between observations and dissipationless simulations, the problem arises because we expect 
that dissipationless simulations predict the same density profile of a system whose physics 
is not just that of dark matter. This lead to the conclusion that the so called cusp/core 
problem is not a reason to reject the CDM model. 



4. Conclusions 



In this paper, we studied the cusp/core problem by means of an improved version of the 
SIM, taking into account, simultaneously and for the first time, the effects of ordered and 
random angular momentum, dynamical friction and adiabatic contraction. Initial conditions 
were introduced by means of the theory of Gaussian random fields. Angular momentum 
was calculated through the standard theory of acquisition of angular momentum through 
tidal torques, while the random part of angular momentum was assigned to protostructures 
according to Avila- Reese et al. (1998) scheme as modified by Ascasibar et al. (2003). Dy- 
namical friction was calculated dividing the gravitational field into an average and a random 
component generated by the clumps constituting hierarchical universes. The adiabatic con- 
traction was taken into account by means of Gnedin et al. (2004) model and Klypin et al. 
(2002) model taking also account of exchange of angular momentum between baryons and 
dark matter. The improved SIM of the present paper, taking account the previous effects 
gives rise to haloes being characterized by log-log density slope that changes gradually within 
the virial radius and slopes of the inner part of haloes flattening with decreasing mass. The 
density proflles of structure having masses smaller than IO^^Mq are well fltted by Burkert's 
profiles. We then calculated the time evolution of a dwarf galaxy having mass IO^Mq. The 
result showed an initial steepening of the profile due to adiabatic compression 
and a subsequent flattening till a slope a ~ is reached, due to angular momen- 
tum and dynamical friction. If the effects of angular momentum and dynamical 
friction are not taken into account, the final profile is cuspy. We then compared 
some of the rotational curves given by Gentile et al. (2004) with the rotational curves ob- 
tained by means of our model, obtaining a good agreement. In the case of clusters of galaxies 
the density profile evolution is similar to that observed on galactic scales with the difference 
that the final slope is steeper than in the dwarf galaxies case. However the total mass profile 
is still cuspy. The behavior of the dark matter halo is in agreement with the analysis of Sand 
et al. (2002, 2004) who fitted the baryonic and dark matter profiles only in the very inner 
part of the cluster MS 2137-23 within ~ 50h'^ kpc by means of a generalized NFW profile. 
The behavior of total mass is in agreement with X-ray observations by Chandra and XMM 
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(Buote 2003, 2004; Lewis et al. 2003) weak lensing (Dahle et al. 2003) and strong lensing 
(Bartelmann 2002), which are consistent with a cusp having a = 1 or larger. The previous 
results tell us that the apparent disagreement between dissipationless N-body simulations 
and observations of the central density profiles of galaxies and clusters can be resolved within 
the CDM paradigm (cosmogony). The cusp/corc problem is simply due to the fact that we 
are trying to predict the dynamics of galaxies at small scales by using just the CDM model 
while at that scales other effects influence their dynamics. 

We would like to thank Massimo Ricotti, Nicos Hiotelis, and Antonaldo Diaferio for 
their very helpful suggestions and comments. 
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A. Model details 

As shown by spherical collapse mode@ of Gunn & Gott (1972), a bound mass shell 
having initial comoving radius Xi will expand to a maximum radius Xm (named apapsis or 
turn-around radius Xta)'- 

Xm = g{Xi) = Xi/6i (Al) 

where the mean fractional density excess inside the shell, as measured at current epoch to, 
assuming linear growth, can be calculated as: 

S^ = -, r 5{y)yHy (A2) 
Jo 

At initial time ti and for a Universe with density parameter fij, a more general form of Eq. 
(lATD (Peebles 1980) is : 

The last equation must be regarded as the main essence of the SIM. It tells us that the final 
time averaged radius of a given Lagrangian shell does scale with its initial radius. Expressing 
the scaling of the final radius, x, with the initial one by relating x to the turn around radius, 
Xm, it is possible to write: 

X = f{Xi)Xm (A4) 

where / depends on a: 

f = f(a) = 0.186 + 0.156a + O.OlSa^ + 0.017a^ - 0.0045a^ + 0.0032a^ (A5) 

(Zaroubi, Naim & Hoffman 1996). If mass is conserved (i.e., m{xi) = m{xm)) and each shell 
is kept at its turn-around radius (/ = 1), the shape of the density profile is given by (Peebles 
1980; HS; White & Zaritsky 1992): 



— g{Xi) — Xi= 7T (A3) 



Pta{Xm) = Pi{Xi) ( j (A6) 

\Xm / (^XfYi 

Using the Virial theorem one obtains a value for the collapse factor / = 0.5 and the final 
density profile is obtained using again mass conservation: 

p{x)x^dx = pixidxi (A7) 



^^A slightly overdense sphere, embedded in the Universe, is a useful non-linear model, as it behaves 
exactly as a closed sub-universe because of Birkhoff 's theorem. The sphere is divided into spherical "shells" . 
A spherical "shell" may be defined as the set of particles at a given radius that are all at the same phase in 
their orbits (see Le Delliou & Henriksen 2003). 
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If one approximates the initial density of the shell of radius Xi as: 

Pi{xi) = pb,i[l + Si{xi)] (A8) 

and expand the right hand-side of Eq. flA3l) in 6i keeping only the linear term, one obtains 
the power-law density profile obtained by HS: 

p{x) oc a;-3("+3)/(-+4) (A9) 



However, after reaching maximum radius, a shell collapses and will start oscillating 
and it will contribute to the inner shells and so even energy is not an integral of motion 
anymore and the collapse factor, /, is no longer constant. The effect of the infalling outer 
shells on the dynamics of a given shell can be described as follows. As I previously told, 
one assumption of SIM is that the collapse is "gentle" . One can assume that the potential 
well near the center varies adiabatically (Gunn 1977, FG84). This means that a shell near 
the center makes many oscillations before the potential changes significantly (Gunn 1977, 
FG84) or similarly the orbital period of the inner shell is much smaller than the collapse time 
of the outer shells (Zaroubi & Hoffman 1993). This implies that the radial action § v{r)dr 
(being v{r) the radial velocity) is an adiabatic invariant of the inner shell. As the outer 
shells collapse, the potential changes slowly and because of the above adiabatic invariant the 
inner shell shrinks. If a shell has an apapsis radius (i.e., apocenter) Xm and initial radius 
Xj, then the mass inside Xm is obtained summing the mass contained in shells with apapsis 
smaller than Xm (permanent component, m^) and the contribution of the outer shells passing 
momentarily through the shell Xm (additional mass madd) ■ Because of mass conservation, we 
have: 

4 _ 

mp{xm) = m{xi) = -TTpbAXi{l + 5i) (AlO) 

where pb^i is the constant density of the homogeneous Universe at the initial time. The 
additional component madd{xm) is: 

dmix) 

rriaddixm) = I PrA^)—^dx (All) 



where R is the radius of the system (the apapsis of the outer shell) and the distribution of 
mass m{x) = m{xm) is given by Eq. ( 1A6I) . The total mass is so given by: 



mT{Xm) = mp{xm) + madd{Xm) (Al2) 

while Px^[x) is the probability to find the shell with apapsis x inside radius x^, calculated 
as the ratio of the time the outer shell (with apapsis x) spends inside radius Xm to its period. 
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This last quantity can be expressed as 

P.. 



■J Xn 



drj 



iv) 



rx 

J Xr, 



drj 



(A13) 



where Xp is the peri center of the shell with apsis x and Vx{ri) is the radial velocity of the shell 
with apapsis x as it passes from radius 77. This radial velocity can be obtained by integrating 
the equation of motion of the shell: 



dVr 



, dr K 



(A14) 



where /i(r, u) is the ordered specific angular momentum generated by tidal torques, j(r, u) 
the random angular momentum (see RG87 and the following of the present paper), G{r) 
the acceleration, A the cosmological constant and /x the coefficient of dynamical friction. I 
shall discuss in the Appendices C and D, how to calculate angular momentum and /i. In the 
peculiar case of /i = 0, Eq. flA14p can be integrated to obtain the square of velocity: 



v{ry 



G 



y"^ Jo y^~" ' 6 



(A15) 



where e is the specific binding energy of the shell that can be obtained from the previous 
equation at turn-around when, dr/dt = 0. 



If /i 7^ 0, the previous equation must be substituted with: 



df 
dt 



TUT A ■ 

y,Z 3 



(A16) 



which can be solved numerically for v. 



Following Gunn (1977) and FG84, the collapse factor /(xj) of a shell with initial radius 
Xi and apapsis Xm is given by: 



m ] 



(A17) 



and the final density profile can be obtained using Eq. ( ]A6I) together with Eq. ( ]A4I) as: 



p{x) 



Pta[,Xr, 
P 



1 + 



d\nf 
ding 



(A18) 



^•^As defined in Appendix B, = d{0)/a, where cr is the mass variance filtered on a scale Rf 
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The collapse factor f{xi) as previously reported and as confirmed by N-body simulations 
(Voglis et al. 1995), is not a constant but it is related to the initial profile of the density 
perturbation and moreover increases with the initial radius. The radial collapse is recovered 
for / — >• and Xj — > (Lokas 2000). The calculation of the collapse factor (Eq. IA17P requires 
that we find before the mass rriadd which namely means to evaluate the integral in Eq. (1A11|) . 



After changing varibles from the turn- around radius to the initial one (Lokas 2000, Hiotelis 
2002), one can calculate it numerically. By means of the quoted change of variables, Eq. 
(lAlip can be written as: 

rriaddirm) = ^T^Pb,i / PxXO[^ + (A19) 

Jxi 

where P^X^'i) = I{xi)/I{x[) with 

/M= r^--^^^d„, (A20) 

= g{xi), Xp = g~^{xp), and where Xp is the pericenter of the shell with initial radius x[. 
The upper limit Xb of the integral is taken to be the initial radius of the sphere that has 
collapsed at the present epoch. The same change of variables must be applied to Eq. ( 1A16I) . 
So, the radial velocity f of a shell with apapsis x = g{xi) as it reaches the radius r = g{ri) 
is given by: 



dvlir) „ o ~\hl+il ^ , , A 



dt 



^ r + -r 
3 



Vx{r) (A21) 



The potential \E' after the change of variables is given by the expression: 

n9(n)\ = ^ + G /" ^^Mflld.,, (A22) 

g{xh) Jr, g^[Xi) dxi 

where the distribution of mass m{xi) is that at the initial conditions. 

Summarizing, knowing the initial conditions, angular momentum and the coefficient of 
dynamical friction, for a given shell one integrates the equation of motions (Eq. IA2ip . then 



one calculates the probability Px^ (Eq. IA20p . from this the contribution of the shell to rriadd 



(Eq. IA19[) . the collapse factor (Eq. IA17I) and finally the density profile through Eq. ( 1A18I) 



(see also Lokas 2000 (Sect. 4); Ascasibar et al. 2004 (Sect. 2.1)). 



B. Initial conditions 



In order to calculate the density profile it is necessary to calculate the initial overdensity 
5i{xi). This can be calculated when the spectrum of perturbations is known. It is widely 



-36- 



accepted that structure formation in the universe is generated through the growth and col- 
lapse of primeval density perturbations originated from quantum fluctuations (Guth & Pi 
1982; Hawking 1982; Starobinsky 1982; BBKS) in an inflationary phase of early Universe. 
The growth in time of small perturbations is due to gravitational instability. The statistics 
of density fluctuations originated in the inflationary era are Gaussian, and can be expressed 
entirely in terms of the power spectrum of the density fluctuations: 



P{k) = (141 



where 



5k — J d^kexp{—i\ac)5{x) 

- Pb 



5(x) 



Pb 



(Bl) 



(B2) 



(B3) 



and Pb is the mean background density. In biased structure formation theory it is assumed 
that cosmic structures of linear scale Rf form around the peaks of the density field, 5(x), 
smoothed on the same scale. HS suggested that, according to the hierarchical scenario of 
structure formation, haloes should collapse around maxima of the smoothed density field 
(see below). The statistics of peaks in a Gaussian random field has been studied in the 
classical paper by BBKS. A well known result is the expression for the radial density profile 
of a fiuctuation centered on a primordial peak of arbitrary height u: 



■m 



-1/2 



(B4) 



e(0)V2 ^(1-^2) 

(BBKS; RG87), where u — S{0)/a (see the following for a definition of a) is the height of a 
density peak, ^ (r) is the two-point correlation function: 

1 



27rV ./o 



P{k)ksm{kr)dk 



7 and R^ are two spectral parameters given respectively by: 

Jk^P{k)dk 



7 



[/ k^P{k)dkJk^P{k)dk] 



1/2 



R* 



3 / k^P{k)dk' 
Jk^P{k)dk 



nl/2 



(B5) 



(B6) 



(B7) 



while 'J?(7i^, 7) is: 



e{u^, 7) = 



3(1 -7') + (1.216- 0.97^) exp -(i)(^) 



3(l-72)+0.45+(^)' 



1/2 



(B8) 
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Then 5i is calculated from Eq. flB4l) similarly to Ascasibar et al. (2003) (their Section 
2.2). In order to calculate 5{r) we need a power spectrum, P{k). The CDM spectrum used 
in this paper is that of BBKS (equation (G3)), with transfer function: 

T{k) = ^ VJ"^^^^^ ■ [1 + 3.89g + (16.1g)2 + (5.46^)^ + (6.71)^'^' (B9) 

where q = j^^fe^Mpc-^ • Here 6 = per/(l-68p^) represents the ratio of the energy density in 
relativistic particles to that in photons {6 = 1 corresponds to photons and three flavors 
of relativistic neutrinos). The spectrum is connected to the transfer function through the 
equation: 

P{k) = PcDMe-"^'"''"f (BIO) 
where -R/ is the smoothing (filtering) scale and Pcdm is given by: 

PcDM = AkT\k) (Bll) 

where A is the normalization constant. We normalized the spectrum by imposing that the 
mass variance of the density field 

1 r°° 

a^M) = — / dkk^P{k)W^{kR) (B12) 
convolved with the top hat window 







W{kR) = -^{smkR-kRcoskR) (B13) 
[kRy 

of radius 8 h^^ Mpc^^ is erg = 0.76 (Romano-Diaz et al. 2008). Throughout the paper we 
adopt a ACDM cosmology with WMAP3 parameters, fi^ = 1 - = 0.24, = 0.76, 
Qb = 0.043 and h = 0.73, where h is the Hubble constant in units of 100 km Mpc^^. 

The mass enclosed in Rf is calculated, as in RG87, as M = An/SpbR}, so that for 
Rf = 0.12 Mpc, M ~ IO^Mq B Structure like Galaxies form from high peaks in the 
density field, high enough so that they stand out above the "noise" and dominate the infall 
dynamics of the surrounding matter. Density profile around particles located at the local 
maxima and minima of the density field is given by < — (x) >=< 6{x) >= z^^(x)/^(0)^/^ 
(Peebles 1984), where ^ is the two-point mass correlation function, and x is the comoving 



^''For precision sake, the mass scale M is connected to the smoothing scale by: Mq = (27r)'^3/2)Ry for a 

Gaussian smoothing {P{k,Rf) = e~^f^ P{k)) and by Mth = 47r/3R|^^ for top hat smoothing. The mass 
enclosed by the smoothing function applied to the uniform background is the same for i?/ = Q.QARth (see 
BBKS). 
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separation. Eliminating minima, the density excess around centers of local density peaks 
was derived by BBKS, and is given by Eq. flB4p . 

The amplitude of any given peak is expressed in terms of its a deviation, where a = 
^(0)^/^. Thus the central density contrast of an ua peak is i^^(O)^''^ and the peak height 
is given hj u = 6{0)/a. Given that galaxies are rather common, they must have formed 
from peaks that are not very rare, say, 2-4 a peaks (RG87)). In Fig. 6, we plot density 
profiles for u = 2,3, and 4. Light lines are density runs around maxima or minima, and 
are proportional to ^{x). Heavy lines show density profiles calculated by BBKS. Note 
that the plotted quantity is the initial excess density distribution around a peak 
linearly evolved to the present day, namely 6{x) = 6iD{zi), where D{z) describes 
the growth of density perturbations (Peebles 1980). As previously reported in Section 
2, the assumption that the initial density profile is proportional to the correlation function 
(light lines in Fig. 6), instead of using the BBKS initial density profile (heavy lines in Fig. 
6) is one of the reasons causing the discrepancies between the SIM and some high resolution 
N-body simulations (DP2000). 



C. Calculation of the angular momentum 
C.l. Ordered angular momentum 

The explanation of galaxies spins gain through tidal torques was pioneered by Hoyle 
(1949). Peebles (1969) performed the first detailed calculation of the acquisition of angular 
momentum in the early stages of protogalactic evolution. More recent analytic computations 
(White 1984, Hoffman 1986, R88; Eisenstein & Loeb 1995; Catelan & Theuns 1996a, b) 
and numerical simulations (Barnes & Efstathiou 1987) have re-investigated the role of tidal 
torques in originating galaxies angular momentum. 

Following Eisenstein & Loeb (1995), we separate the universe into two disjoint parts: 
the collapsing region, characterized by having high density, and the rest of the universe. The 
boundary between these two regions is taken to be a sphere centered on the origin. As usual, 
in the following, we denote with p(x), being x the position vector, the density as function 
of space and 5(x) = , The gravitational force exerted on the spherical central region 

by the external universe can be calculated by expanding the potential, $(x), in spherical 
harmonics. Assuming that the sphere has radius R, we have: 

°° Ayr ' 

'^(^) = E ^ E «lm(^)>1m(^, 0)X' (CI) 
/=0 m=-l 
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where Y|m are spherical harmonics and the tidal moments, aim, are given by: 

POO 

aUx)=pb rim(e,0)p(s)s-'^^rf=^s (C2) 
Jr 

In this approach the proto-structure is divided into a series of mass shells and the torque 
on each mass shell is computed separately. The density profile of each proto-structure is ap- 
proximated by the superposition of a spherical profile, S{r), and a random CDM distribution, 
e(r), which provides the quadrupole moment of the proto-structure. To the first order, the 
initial density can be represented by: 

p(r)=pb[l + 5(r)][l+e(r)] (C3) 

where ^(r) is given by: 

{\ek\') = P{k) (C4) 

being P{k) the power spectrum. The torque on a thin spherical shell of internal radius x is 
given by: 

CM r 

t{x) = -—-^ / £(x)xxv$(x)dfi (C5) 



47r 

where Mgh = 47rpb [1 + ^{x)] x'^6x. Before going on, 1 want to recall that we are interested 
in the acquisition of angular momentum from the inner region, and for this purpose we take 
account only of the / = 2 (quadrupole) term. In fact, the / = term produces no force, while 
the dipole (/ = 1) cannot change the shape or induce any rotation of the inner region. As 
shown by Eisenstein & Loeb (1995), in the standard CDM scenario the dipole is generated 
at large scales, so the object we are studying and its neighborhood move as bulk flow with 
the consequence that the angular distribution of matter will be very small, then the dipole 
terms can be ignored. Because of the isotropy of the random field, e(x). Equation (ICSp can 
be written as: 

< kr >= \/(30)^ [< a2mixf >< q2m{xf > - < a2™(x)g2*„(x) >^Y^^ (C6) 

where <> indicates a mean value of the physical quantity considered. 

In order to find the total angular momentum imparted to a mass shell by tidal torques, it 
is necessary to know the time dependence of the torque. This can be done connecting ^2111 and 
a2m to parameters of the spherical collapse model (Eisenstein & Loeb 1995 (equation (32), 
R88 (equation (32) and (34)). Following R88 we have: 

^ ^-3 (l-cosg)V2(g) 

q2m{0) = -;q2mfiOo T"^ (C7) 

/lW-(|)/2W 
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and 



4X 4/3_ 



-4 



a2m(^) = a2n.,o ( 3 1 5o{9 - sine)- (C8) 



The collapse parameter 6 is given by: 

t{e) = ^to6o''^\e-sme) (C9) 

Equation (1C7I1 and (ICSj) . by means of equation (1C6I) . give to us the tidal torque: 

where fi(9) and /2(6') are given in R88 (Eq. 31), tq and 60 = are respectively the torque 
and the mean fractional density excess inside the shell, as measured at current epoch to- The 
angular momentum acquired during expansion can then be obtained integrating the torque 
over time: 

As remarked in the Del Popolo et al. (2001) the angular momentum obtained from Eq. 
flC12p is evaluated at the time of maximum expansion tm. Then the calculation of the 



angular momentum can be solved by means of Eq. flC12p . once we have made a choice 
for the power spectrum. With the power spectrum (filtered on a galactic scale) and the 
parameters given in Appendix A, for a. u = 2 peak of mass ~ 2 x lO^^M©, the model gives 
a value of 2.5 x 10^^ 

s 

It is interesting to compare with a different method like that of Catelan & Theuns 
(1996), who calculated the angular momentum at maximum expansion time (see their Eqs. 
(31)- (32)) and compared it with previous theoretical and observational estimates. Assuming 
the same value of mass u used to obtain our previously quoted result same distribution of 
angular momentum (//) as adopted by Catelan & Theuns (1996), we obtain a value for the 
angular momentum of 2.0 x 10''^ ^^j- in agreement with our result. 

For what concerns the calculation of the random angular momentum, j, it was described 
in next subsection. 



C.2. Random and ordered angular momentum 



As reported in the introduction, several authors have emphasized the effect of an 
isotropic velocity dispersion (thus of non-radial motion) in the core of collisionless haloes 
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(RG87; White & Zaritski 1992; Avila-Reese et al. 1998; Hiotelis 2002; Nusser 2001; Le 
Delliou & Enriksen 2003; Ascasibar et al. 2004; Williams et al. 2004). One common result 
of the previous studies is that larger amount of angular momentum leads to shallower final 
density profiles in the inner region of the halo. There are two sources of angular momentum 
of coUisionless dark matter: (a) bulk streaming motions, and (b) random tangential motions. 
The first one (ordered angular momentum (RG87)), as described in the previous subsection, 
arises due to tidal torques experienced by proto-halos, and is usually quantified as a di- 
mensionless spin parameter A (Peebles 1969). The second one (random angular momentum 
(RG87)) is connected to random velocities (see RG87 and the following of the present paper). 

In the present paper, we took into account both types of angular momentum: random 
j, and ordered, h. Type (a), connected to tidal torques is obtained as described in Appendix 
CI, obtaining the rms torque, T(r), on a mass shell using Eq. ( ]C5[) and then calculating 
the total specific angular momentum, h{r, u), acquired during expansion by integrating the 
torque over time (Ryden 1988a (hereafter R88), Eq. 35): 

where Mgh = 4'7rpb [1 + ^i^)] x^5x is the mass in a thin spherical shell of internal radius x. 
To is the tidal torque at time to, the functions fi{^), fii,'^) are given by R88 (Eq. 31) while 
the mean over-density inside the shell, 5(r), is given by Eq. ( ]A2[) . 

In Fig. 7, we show the variation of /i(r, v) with mass for three values of the peak 
height z/ = 2, 3, 4 for a power spectrum filtered on scale Rf = 0.12 Mpc. The rms specific 
angular momentum, h{r,v), increases with distance r while peaks of greater v acquire less 
angular momentum via tidal torques. This is the angular momentum-density anti-correlation 
showed by Hoffman (1986). This effect arises because the angular momentum is proportional 
to the gain at turn around time, t^, which in turn is proportional to (5(r, i/)~^/^ oc z/^^/^. 

For what concerns the calculation of the random angular momentum, j, we have to recall 
that this quantity is product of highly non-linear processes and it is difficult to calculate. An 
estimate has been performed by RG87 and R88 assuming that the substructure of a peak can 
be described as the sum of the mean spherical profile, 5(x), and of a random Gaussian field 
with the same power spectrum as the parent density field (see Appendix C, Eq. IC3p . RG87 
showed that j(r) = (2/3)^/^rmAf (r, tm/2) where Af is the r.m.s. velocity. It is important 
to note that j, as /i, are anti-correlated with the height of the peak since both Af and r^ax 
decrease with an increase in the height v of the peak. To the lowest order in r/Rf the 



^■^As previously reported radius scales with mass as r oc M^f^ . 
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random angular momentum j is (Ryden 1988b): 

(C13) 



j{r) 1.5 r 2 



The ratio of ordered to random angular momentum, to the same order is: 

h 0.09 , r .0, 

The random angular momentum is important for the particle (shell) orbits since the larger 
it is, the larger is the orbital ellipticity, and then the shell penetrates less to the center, 
resulting in a flattening of the inner density profile. In several of the previously quoted 
papers only the second type of angular momentum (the one generated by random motions) 
was taken into consideration. The usual approach (Nusser 2001; Hiotelis 2002; Ascasibar et 
al. 2003) consists in assigning a specific angular momentum at turn around: 

3 oc \/GM^ (C15) 

With this prescription, the orbital eccentricity e is the same for all particles in the halo 
(Nusser 2001). Avila- Reese et al. (1998) expressed the specific angular momentum j through 
the ratio Cq = ( ^^^^ ) , and left this quantity as a free parameter. Processes related to 
mergers and tidal forces that could produce tangential perturbations in the collapsing matter 
were implicitly considered in their model trough the free parameter cq. According to them, 
the detailed description of these processes is largely erased by the virialization process, 
remaining only through the value of cq, which then they fixed to cq = 0.3. The value 
e ~ 0.2 gives density profiles very close to the NFW profile (Avila- Reese et al. 1998, 1999). 
This procedure is justified by N-body simulations of halo collapse: for CDM halos, N-body 
simulations produce constant < ^-^^ >~ 0.2 ratios of dark matter particles in virialized 
, and the range for this value is around 0.1-0.3 (e.g., Ghigna et al. 1999). Ascasibar 
et al. (2003) noticed that particle orbits are slightly more radial as we move out to the current 
turn around radius, rta- There is a dependence on the dynamical state: major mergers are 
well described by constant eccentricity up to the virial radius 0, while relaxed systems are 
more consistent with a power-law profile. Minor mergers are in the middle. The average 
profile can be fitted by a power law, but the slope is shallower than for relaxed systems. A 
least-square fit to the relaxed population yields: 

e(w) ~ 0.8(w/n„)°-i (C16) 



^^The constancy of < > throughout the halo has been interpreted as a proof that the adiabatic 

approximation is largely valid for haloes generated in N-body simulations (Jesscit et al. 2002). 

^^For their dark matter haloes, r^/rta is typically of the order of 0.2-0.3. 
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for Tmax < O.lrfQ. In the present paper, we use Avila-Reese et al. (1998) method with the 
correction of Ascasibar et al. (2003). It is important to notice that even using a constant 
value for the eccentricity angular momentum will change with radius and v. In fact the 
specific angular-eccentricity relation is: 

J = ^GMr^^l-e) (C17) 

which depends from mass (radius) and is anti-correlated with v (r^ oc 

In Fig. 8, we plot the specific random angular momentum for different values of v 
showing a similar behavior to that of the ordered angular momentum. Since the ordered 
angular momentum h oi a, shell increases more rapidly with the initial radius of the shell 
than random angular momentum j (R88), the inner regions of the halo are "hotter" than 
the outer regions, in the sense of having a smaller ratio of ordered to random velocities. 
Moreover high peaks are hotter than low peaks. This behaviour is plotted in Fig. 9, where 
we plot the ratio of ordered to random angular momentum. 

It is important to note that after turn-around, not only the random angular momentum 
i will contribute to originate non-radial motions in the protostructure but also the ordered 
angular momentum. In fact, as successive shells expand to their maximum radius they 
acquire angular momentum and then fall in on orbits determined by the angular momentum. 
Actual peaks are not spherical; thus the infall of matter will not be purely radial. Random 
substructure in the region surrounding the peak will divert infalling matter onto non-radial 
orbits. The role of this random motions is of fundamental importance in the structure 
formation. 

In the present paper, as in Nusser (2001), and Williams et al. (2004)0, the contribution 
of h and j are evaluated analytically but are not imparted to the shell until it reaches turn- 
around. In other words, as long as a given shell is expanding its dark matter particles' 
positions and velocities are not corrected for the effects of h and j; it is assumed that 
during this relatively orderly evolutionary stage the effects of random perturbations are 
small compared to what they will be later when the halo starts to collapse. Finally, angular 
momentum is taken into account by adding the term ^ in the equation of acceleration of 
the shell. 



'Note that in those papers, just the random angular momentum is taken into account. 
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D. Dynamical friction. 



The effect of dynamical friction on the central halo profile has been studied in several 
papers (e.g., EZOl; Tonini, Lapi & Salucci 2006 (TLS)) using different approaches. EZOl 
studied the changes induced by dynamical friction in the dark matter halo structure within 
the context of galaxy formation, where gas cools to form dense clumps. The orbital energy 
lost by the clumps to the dark matter background is sufficient to "heat up" and flatten 
the dark matter density cusps. TLS performed two independent calculations over the NFW 
halo: 1) they perturbed its phase-space distribution function with angular momentum and 
2) they followed the evolution of the baryonic angular momentum during the collapse inside 
the halo potential well. In both cases the final result was a core-like profile. In the present 
paper, we took into account dynamical friction by introducing the dynamical friction force 
(see the following) in the equation of motions. 

In a hierarchical structure formation model, the large scale cosmic environment can 
be represented as a collisionless medium made of a hierarchy of density fiuctuations whose 
mass, M, is given by the mass function N{M,z), where z is the redshift. In these models 
matter is concentrated in lumps, and the lumps into groups and so on. In such a material 
system, gravitational field can be decomposed into an average field, Fo(r), generated from 
the smoothed out distribution of mass, and a stochastic component, FstochiT^), generated 
from the fluctuations in number of the field particles. The stochastic component of the 
gravitational field is specified assigning a probability density, W(F), (Chandrasekhar & von 
Neumann 1942). In an infinite homogeneous unclustered system H^(F) is given by Holtsmark 
distribution (Chandrasekhar & von Neumann 1942) while in inhomogeneous and clustered 
systems W^(F) is given by Kandrup (1980) and Antonuccio-Delogu & Barandela (1992) 
respectively. The stochastic force, Fstoch, in a self-gravitating system modifies the motion 
of particles as it is done by a frictional force. In fact a particle moving faster than its 
neighbors produces a deflection of their orbits in such a way that average density is greater 
in the direction opposite to that of traveling causing a slowing down in its motion. Following 
Chandrasekhar & von Neumann's (1942) method, the frictional force which is experienced 
by a body of mass M (galaxy), moving through a homogeneous and isotropic distribution of 
lighter particles of mass m (substructure), having a velocity distribution n{v) is given by: 



where log A is the Coulomb logarithm, p the density of the field particles (substructure). 
A more general formula is that given by Kandrup(1980) in the hypothesis that there are no 
correlations among random force and their derivatives: 



M-^ - -A7rG^M^n{v)^logAp 




(Dl) 



F = - 



/ W{F)F^T{F)d^F 
2<v^> 



V 



(D2) 
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where F is dynamical friction force per unit mass, fi is the coefficient of dynamical friction, 
T{F) the average duration of a random force impulse, < f ^ > the characteristic speed of a 
field particle having a distance r ~ from a test particle (galaxy). This formula is 

more general than Eq. ( 1D1[) because the frictional force can be calculated also for inhomo- 
geneous systems when W{F) is given. If the field particles are distributed homogeneously 
the dynamical friction force per unit mass is given by: 

(Kandrup 1980), where ma and are respectively the average mass and number density of 
the field particles. Using the virial theorem we also have: 



where Mt^t is the total mass of the system, Rgys its radius and is the total number of field 
particles. The dynamical friction force per unit mass can be written as follows: 

4.44[Gman„J^/2 . v v , , 

F = -f,v = L_p^ log {i.i2iV^/3} = -e,-^ (D5) 

where = ^R^^y^na, a is the expansion parameter, connected to the proper radius of a shell 
by: 

r(ri,t) = ria{ri,t) (D6) 

nac = Ua X is the comoving number density of field particles. This last equation supposes 
that the field particles generating the stochastic field are virialized. This is justified by the 
previrialization hypothesis (Davis & Peebles 1977). To calculate the dynamical evolution of 
the galactic component of the cluster it is necessary to calculate the number and average 
mass of the field particles generating the stochastic field. 

The protocluster, before the ultimate collapse at z ~ 0.02, is made of substructure having 
masses ranging from 10^ — IO^Mq and from galaxies. I suppose that the stochastic gravita- 
tional field is generated from that portion of substructure having a central height u larger 
than a critical threshold z/c. This latter quantity can be calculated (following Antonuccio- 
Delogu & Colafrancesco 1994 (hereafter ADC)) using the condition that the peak radius, 
^pki^ > i^c), is much less than the average peak separation na{i^ > Uc)^^^^, where is given 
by the formula of BBKS for the upcrossing points: 



exp(z/2/2) ^ 7 



n.(.>.J---^(^)^k^-l + 

■exp(-57V,Vl8)] (D7) 



572(1-572/9)1/2 
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where 7, i?* are parameters related to moments of the power spectrum (BBKS Eq. 4.6A). 
The condition rpfc(z/ > u^,) < O-lUaiu > i^c)~^^^ ensures that the peaks of substructure are 
point hke. Using the radius for a peak: 



V2R, 



^^fTo)(73 + (0.9/z/))3/2 



1/3 



(D8) 



(ADC), I obtain a value of z/^ = 1.3 and then we have na(z/ > Uc) = 50.7 Mpc (7 = 0.4, 
-R=K = 50 kpc) and is given by: 

1 f°° 

ma = , . / mpk{v)Npk{v)dv = lO^Mo (D9) 

(in accordance with the result of ADC), where rupk is given in Peacock & Heavens (1990) 
and Np}. is the average number density of peak (BBKS Eq. 4.4). Galaxies and Clusters of 
galaxies are correlated systems whose autocorrelation function, ^(r), can be expressed, in a 
power law form (Peebles 1980; Bahcal & Soneira 1983; Postman et al. 1986; Davis & Peebles 
1983; Gonzalez et al. 2002). The description of dynamical friction in these systems need to 
use a distribution of the stochastic forces, W{F), taking account of correlations. In this last 
case the coefficient of dynamical friction, /x, may be calculated using the equation: 

^ = j d^YW{F)F^T{F)/{2 < >) (DIO) 

and using Antonuccio-Delogu & Atrio-Barandela (1992) distribution: 

1 f°° 

W{F) = / dkksin{kF)Af{k) (Dll) 

where Af, which is a linear integral function of the correlation function ^(r), is given in the 
quoted paper (Eq. 36). 

In Fig. 10, we plot the coefficient of dynamical friction of a system having Rsys = 5h~^ 
Mpc for different values of and Rj. The increase of /i with is mainly due to an 
increase of the average mass of the peaks generating the stochastic field and a simultaneous 
decrease of their number density. 



In Fig. 11, we show the effects of non-radial motions (angular momentum) and dynam- 
ical friction, in the case of a z/ = 3 peak. We plot the evolution of the expansion parameter, 
a(t), with time expressed in terms of the collapse time of a pure SIM model (namely when 



^^We calculate the contribution to /i from small peaks of the density field, namely those peaks small 
enough to be considered "pointlike", having their central height v larger than a critical threshold i/c- 



tidal torques, dynamical friction and cosmological constant are not taken into account). As 
shown by Gunn & Gott (1972), this last quantity is given by: 

Troir, u) = — = 

^ ' H[5{r,y)fl^ 

where H is Hubble constant. 

In Fig. 11, the dashed line represents the effect of dynamical friction while the dotted 
line the cumulative effect of dynamical friction and angular momentum. The presence of 
non-radial motions and dynamical friction changes the dynamics of the shells producing an 
increase in the time of collapse. The collapse delay is larger for low value of v and reduces 
for high values of v for which the collapse time becomes equal to that of the pure spherical 
collapse. This result is in agreement with the angular momentum-density anticorrelation 
effect: density peaks having low value of u acquire a larger angular momentum than high u 
peaks and consequently the collapse is more delayed with respect to high v peaks. Dynamical 
friction has a similar effect to that of angular momentum for shell evolution. 
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E. Baryonic dissipative collapse 

The response of dark matter to baryonic infall has traditionally been calculated using 
the model of adiabatic contraction. Eggen et al. (1962) were the first to use adiabatic 
invariants of particle orbits to estimate the effect of a changing potential in a contracting 
proto-galaxy. Zeldovich et al. (1980) presented the first analytical expressions for adiabatic 
contraction (AC) for purely radial and circular orbits, as well as the numerical tests of such 
model. The present standard form of the AC model was introduced and tested numerically 
by Blumenthal et al. (1986) (see also RG87). Blumenthal et al. (1986) described an 
approximate analytical model which has been checked by numerical simulations (Barnes 
1987; Oh 1990) for calculating the radial redistribution of the dissipationless halo matter of 
a protostructure when its dissipational matter falls in toward the center. One starts by noting 
that for a particle moving in a periodic orbit § pdq (where p is the canonical momentum 
conjugate to the coordinate q) is an adiabatic invariant. Provided that M(r), the mass inside 
the orbital radius r, changes slowly compared with an orbital time (Steigman et al. 1978; 
Zeldovich et al. 1980; RG87) and that particles move in circular orbits about a spherically 
symmetric mass distribution, the mass invariant is rM(r). For purely radial orbits r^M(r^) 
is also an adiabatic invariant, provided that M(r) varies in a self-similar fashion. Consider 
a spherically symmetric protostructure that consists of a fraction F << 1 (this fraction is 
fixed using WMAP collaboration data (Spergel et al. 2003)) of dissipational baryons and 
a fraction 1 — F of dissipationless dark matter particles constituting the halo. One then 
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assumes that the dissipational baryons and the halo particles are well mixed initially (i.e., 
the ratio of their densities is F through the protostructure) . Since there is more phase 
space for nearly circular orbits than for nearly radial orbits, one can assume that the dark 
matter particles move in circular orbits about the proto-structure center with almost random 
oriented angular momentum vectors. As the baryons dissipatively cool and fall into a final 
mass distribution Mb(r), which is constrained by the initial angular momentum distribution, 
a dark matter particle initially at radius rj will move in to radius r < r^. The adiabatic 
invariant for such a particle orbits implies that: 

r [Mb(r) + Mrf„(r)] = nMi^n) (El) 

where Mi{ri) is the initial total mass distribution, Mb{r) (as previously reported) is the final 
mass distribution of dissipational baryons and M^m is the final distribution of dissipationless 
halo particles. Assuming that the orbits of the halo particles do not cross, then 

MdUr) = (1 - F)M,{r,) (E2) 

Eqs. (lEip . (lE2p can be iteratively solved to calculate the final radial distribution of the halo 
particles once Mj(rj) and Mb(r) are given. If F << 1, then for a halo particle not too near to 
the center of the protostructure, the mass interior to its orbit will undergo a small fractional 
change in one orbital period, even if dissipation occurs rapidly provided the particle is fairly 
far from the protostructure center. Therefore, the adiabatic invariant given in Eq. flEip is 
expected to be a good approximation for all but the innermost particles. 

I recall that the AC model has been studied further by R88, Ryden (1991), and Flores 
et al. (1993). It is routinely used in mass modeling of galaxies (Flores et al. 1993, Dalcanton 
et al. 1997, Mo et al. 1998, Courteau & Rix 1999, van den Bosch 2001, van den Bosch & 
Swaters 2001, Klypin et al. 2002, Seljak 2002) and clusters of galaxies (Treu & Koopmans 
2002). 

In order to use the adiabatic invariant, one has to know the initial mass distribution 
Mj(rj) and the final distribution of dissipational baryons Mb(r). A usual assumption is that 
initially baryons had the same density profile as the dark matter (Mo et al. 1998; Cardone 
& Sereno 2005; Treu & Koopmans 2002; Keeton 2001). 

The final baryons distribution is assumed to be a disk (for spiral galaxies) (Blumenthal 
et al. 1986; Flores et al. 1993; Mo et al. 1998; Klypin et al. 2002; Cardone & Sereno 2005). 
In our calculations, we shall assume Klypin et al (2002) model for baryons distribution (see 
their subsection 2.1), when dealing with mass scales typical of spiral galaxies. In the case 
of elliptical galaxies and clusters a typical assumption (that we shall use in the paper) is 
that baryons collapse to a Hernquist configuration (Rix et al. 1997; Keeton 2001; Treu & 
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Koopmans 2002), having density profile: 

P(^) = 1—1 — TTT"^ 1 — ■ (E3) 

(r/r,)(l + r/r,)3 

This profile is described by a scale radius r^, but the projected profiles of elliptical galaxies are 
usually described by an effective (or half-mass) radius Re\ they are related by = 0.551-Re- 
The total mass of a Hernquist model is M = l-npsT^. As previously reported, in general, 
adiabatic contraction is computed numerically, but with a Hernquist model (Eq. IE3P the 
problem can be solved analytically (Keeton 2001). Each initial radius rj maps to a unique 
final radius r given by the solution of the equation 



which is a cubic polynomial in r. Here the galaxy scale radius Vg from Eq. (]E3P has been 
relabeled as Sg. Also, vniir^ = Mi{ri) /M^ is the initial mass profile normalized by the virial 
mass. In the limit r >> s^, Eq. (]E4p has the simple asymptotic solution 

„ - ^^"^^(^^) (E5) 



/ + (l-/)m,(r,) 

The full general solution can also be found analytically, although it cannot be written quite 
so compactly (see Keeton 2001). Once the cubic equation has been solved to map to r, 
Eqs. flEll IE2p can be used to write the total mass profile as 

Mtot{r) = M,{r) + Mdm{r) = ^ M,(rO . (E6) 

This solution of adiabatic contraction by a Hernquist galaxy can be used for any form of the 
initial halo, by simply inserting the desired initial profile Mj(rj) into Eq. (lE4p and Eq. (lE6p . 

Recently Gnedin et al. (2004) presented the first test of the AC model in the self- 
consistent simulations of hierarchical structure formation and propose a simple modification 
which describes numerical results more accurately. They proposed a modified adiabatic 
contraction model based on conservation of the product of the current radius and the mass 
enclosed within the orbit-averaged radius: 

M(f)r = const. (E7) 

where the orbit-averaged radius is 

2 dr 

f= - / r-, (E8) 

-'■r Jrp '^r 

where is the radial period, is the apocenter radius and Vp the pericenter radius. 
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The previous, classical AC model, assumes no angular momentum exchange between 
different components (e.g., baryons and dark matter). Models with exchange of angular 
momentum between the baryons and dark matter are more complicated. The exchange 
probably happens at late stages of the baryonic infall when the baryon density becomes 
large and a non-axisymmetric component may develop due to the excitation of spiral waves 
and/or bar-like modes. Dynamical friction can then result in a transfer of angular momentum 
from the baryons to the dark matter. The classical approach can then be used to compute 
the state of the system at the end of the first stage of adiabatic compression (in which angular 
momentum is conserved). For the later evolution, one can use a simple model of Klypin et 
al. (2002). If one considers a spherical shell of dark matter with radius r, thickness dr, 
density pdm, and specific angular momentum 

J = rV, = [Mb(r) + Md^(r)]r. (E9) 

one can get an implicit equation for the final radius rj: 

(ElO) 



3f = 2 



^ AAM 



AM = Mb,f-Mb. (E12) 

Klypin et al. (2002), where Mb^f is the final baryons mass and M = + Mb is the 
total mass inside a radius r. Eq. flElOp is solved numerically. The solution also gives 
the mass inside a final radius rj. Eq. flElOp has the same structure as Eq. flEip . The 
only difference is the term on the right-hand-side, which is the correction due to angular 
momentum deposition. 
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Fig. 1. — Dark matter haloes generated with the model of Section 2. In panels (a)-(d) 
the solid line represents the NFW model while the dotted line the density profile obtained 
with the model of the present paper for masses 10^ Mq (panel a), lO^^M© (panel b), lO^^M© 
(panel c) and lO^^M© (panel d). The NFW profile for the given mass was calculated with 
Eqs. ([I]),([2]), ([3]) expressing the scaling radius in terms of the virial radius. The dashed 
line in panel (a) represents the Burkert fit to IO^Mq halo. The dashed line in panel (d) 
represents the density profile obtained reducing the magnitude of h, j and fj, as described in 
the text. 
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Fig. 2. — Distribution of the total specific angular momentum, Jtoi- The dotted-dashed and 
dashed line represents the quoted distribution for the halo n. 170 and n. 081, respectively, 
of van den Bosch et al. (2002). The dashed histogram is the distribution obtained from our 
model for the lO^^M© halo and the solid one the angular momentum distribution for the 
density profile reproducing the NFW halo, as descried in Section 4. 
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Fig. 3. — Density profile evolution of a IO^Mq halo. The solid line represents the profile at 
z — 10. The profile at2; = 3,2; = 2,2; = l,2; = 0is represented by the long-dashed line, 
short-dashed hne, dot-dashed hne, dotted hne, respectively. 
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Fig. 4. — Comparison of the rotation curves obtained with the model in Section 2 (sohd 
hnes) with the rotation curves of four LSB galaxies studied by Gentile et al. (2004). The 
dotted line represents the fit with NFW model (see Section 4 for details). 
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Fig. 5. — The density profile evolution of a halo. The (uppermost) dot-dashed 

line represents the total density profile of a IO^^Mq halo at 2; = 0. The profile at 2; = 3, 
z — 1.5, z — 1 and 2; = is represented by the solid line, dotted- line, short-dashed- line, 
long-dashed-line, respectively. 
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Fig. 6. — Initial density profiles of dark matter halos linearly evolved to the present day 
(see text). Heavy lines are density profiles around maxima only (BBKS). Light lines are 
proportional to the two-point correlation function, and represent density run averaged around 
local density maxima or minima in the initial Gaussian random field. Solid, short-dashed, 
and long-dashed lines represent 2, 3, and 4 a peaks. 
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Fig. 7. — The specific angular momentum for tliree values of the parameter i/ (i/ = 2 solid 
hne, u — 3 dotted line, u — A dashed hne) and for Rf — 0.12/i~^Mpc. The radius r is 
connected to the mass M as described in the text. 
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Fig. 8. — The random specific angular momentum, j for tliree values of the parameter u 
{u — 2 solid line, u — 3 dotted hne, u — A dot-dashed line) and for Rf — 0.12h~^Mpc. The 
radius r is connected to the mass M as described in the text. 
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Fig. 9. — The ratio of the ordered, h, to random secific angular momentum, j for three 
values of the parameter u {u — 2 solid hne, u — 3 dotted hne, u — A dot-dashed line) and 
for Rf — 0.12h~^Mpc. The radius r is connected to the mass M as described in the text. 
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Fig. 10. — Variation of the dynamical friction coefficient /j, with i/^ for a system having 
Rsys = 5h^^ Mpc and for Rf = 0.12 Mpc (sohd hne), Rj = 0.3 Mpc (dot-dashed fine), 
Rf = 0.5 Mpc (short-dashed hne), Rf = 0.7 Mpc (long-dashed line). In Eq. (??) units have 
been expressed in terms of Mq, Mpc, T^q. 
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Fig. 11. — The time evolution of the expansion parameter. The sohd hne is a(t) for the pure 
radial collapse (Gunn & Gott 1972); the dashed line is a(t) taking account only dynamical 
friction; the dotted line is a(t) taking account of the cumulative effect of non-radial motions 
and dynamical friction in the case of a = 3 peak. 
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